published in Phys. Rev. E 82, 026706 (2010) 

Stochastic simulation algorithm for the quantum linear Boltzmann equation 

Marc Busse/ Piotr Pietrulewicz,^ Heinz-Peter Breuer,^ and Klaus Hornberger^^'^ 

^Arnold Sommerfeld Center for Theoretical Physics, 

Ludwig-Maximilians-Universitdt Miinchen, Theresienstrafie 37, 80333 Munich, Germany 

'"^ Physikalisches Institut, Universitdt Freiburg, Hermann-Herder-Strafie 3, 79104 Freiburg, Germany 

Max Planck Institute for the Physics of Gomplex Systems, Nothnitzer Strafle 38, 01187 Dresden, Germany 

(Dated: May 3, 2010) 

We develop a Monte Carlo wave function algorithm for the quantum linear Boltzmann equa- 
tion, a Markovian master equation describing the quantum motion of a test particle interacting 
with the particles of an environmental background gas. The algorithm leads to a numerically effi- 
cient stochastic simulation procedure for the most general form of this integro-differential equation, 
which involves a five-dimensional integral over microscopically defined scattering amplitudes that 
account for the gas interactions in a non-perturbative fashion. The simulation technique is used to 
assess various limiting forms of the quantum linear Boltzmann equation, such as the limits of pure 
coUisional decoherence and quantum Brownian motion, the Born approximation and the classical 
Oh' limit. Moreover, we extend the method to allow for the simulation of the dissipative and decohering 

y-» I dynamics of superpositions of spatially localized wave packets, which enables the study of many 

^AJ , physically relevant quantum phenomena, occurring e.g. in the interferometry of massive particles. 
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I. INTRODUCTION 



pH , The motion of a quantum particle interacting with a surroundings particle gas is characterized by collision-induced 

C^ ' decoherence as well as dissipation and thermalization effects. An appropriate master equation which provides a unified 

^ I quantitative description of both phenomena in a mathematically consistent way is the quantum linear Boltzmann 

O^. equation (QLBE), proposed in its weak-coupling form in [l|, |2|, and in final form in [3|. This equation represents 

the quantum mechanical generalization of the classical linear Boltzmann equation which describes the motion of a 

^s) . distinguished test particle under the influence of elastic collisions with an ideal, stationary background gas. The 

^ ' QLBE may be derived on the basis of a monitoringapproach [j] which permits a non-perturbative treatment of the 

l/~) , interactions with the environmental gas particles (3|, |5|- These interactions may therefore be strong and the test 

^^ ' particle may be in a state which is far from equilibrium. A condition for the applicability of the monitoring approach 

^IJ I is that three-particle collisions are sufficiently unlikely, and that successive collisions of the test particle with the 

^^ ■ same gas particle are negligible on the the relevant time scale. These conditions are fulfilled in the case of an ideal 

\^ [ background gas in a stationary equilibrium state. A further condition is that the interactions are short-ranged so that 

scattering theory may be applied. 

The mathematical structure of the QLBE is rather involved and analytical solutions of this equation are known only 

for some specific limiting cases [6|- Moreover, the spatially nonlocal structure of this equation makes a direct numerical 

integration through deterministic methods extremely demanding. However, being in Lindblad form, the QLBE allows 

one to apply the Monte Carlo wave function techniques |7Hl]| . As has been demonstrated in Ref . [l^ these techniques 

; ^ 1 lead to a simple and numerically efficient stochastic simulation method in the momentum representation of the test 

Ci ] particle's density matrix, employing the translational covariance of the QLBE. 

The simulation technique developed in [l2| is restricted to the QLBE within the Born approximation [l|, y|, in 
which the scattering cross section depends only on the momentum transfer of the scattered particles, yielding a 
considerably simplified equation of motion. Here we generalize this stochastic approach to the full QLBE allowing 
for an arbitrary form of the microscopic interaction between the test particle and the ambient gas particles and, 
thus, arbitrary scattering amplitudes. In addition, we show how the algorithm can be extended to simulate efficiently 
the dynamics of spatially localized wave packets. This enables the exact numerical treatment of many physically 
relevant phenomena, such as the loss of coherence in position space and the determination of the fringe visibility in 
interferometric devices, as well as the assessment of the quality of various approximations of the QLBE. 

The paper is organized as follows. Section |TT] contains a brief account of the QLBE and summarizes the most 
important limiting forms of this equation. In Sec. IIIII we develop the Monte Carlo simulation algorithm for the full 
three-dimensional QLBE in momentum space. Our numerical simulation results are presented in Sec. IIVI We discuss 
examples for the decoherence of superpositions of momentum eigenstates, the loss of coherence of superpositions of 
spatially localized wave packets, the decohering influence of the background gas on the fringe visibility of interference 
experiments, relaxation and thermalization processes, and the diffusion limit. Finally, SecFVlcontains a brief summary 
of the results and our conclusions. 



II. THE QUANTUM LINEAR BOLTZMANN EQUATION 

A. General form of the master equation 

The quantum linear Boltzmann equation (QLBE) is a Markovian master equation for the reduced density operator 
p describing the evolution of a test particle in an ideal gas environment. It has the form p — Cp, where the generator 
is of Lindblad structure, 



Cp ^ — 

in 
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^ + HniP),P 



Vp. (1) 



Here iJ„ (P) describes the energy shift due to the interaction with the background gas; it will be neglected in the 
following, since it is usually small. The incoherent part of the interaction is accounted for by the supcropcrator V, 
which can be expressed as 0, |5|, [g 

-Dp = f dQ f dfcx(e*'^''/''i(fc±,P,Q)pi^(fei,P,Q)e-'^-^/'^ 

-i{p,it(fc^,p,Q)L(fc^,P,Q)}), (2) 

with X = (Xi,X2,X3) the position and P = (Pi, P2, P3) the momentum operator of the test particle. The integra- 
tion variables are given by Q, the momentum transfer experienced in a single collision, and k±, corresponding to 
the momentum of a gas particle. The fc_L-integration is carried out over the plane Q = |fc^ G ]R'^|fcx • Q = 0} 
perpendicular to the momentum transfer Q. 

The operator- valued function L {k±, P, Q) contains all the details of the coUisional interaction with the gas; these 
are the gas density ngas, the momentum distribution function p{p) of the gas, and the elastic scattering amplitude 
f{pf,p,). It is defined by da i 



-^TFrf Prcl fc_L,P±Q - —,Prc\{k±,P±Q) + — 
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Here m^, = mM/ {m + M) is the reduced mass, Q = \Q\ gives the modulus of the momentum transfer Q, and the 
function 

Prol P,P = P T7^' 4 

m Ai 

defines relative momenta. The subscripts \\Q and _L Q denote the parts of a given vector P parallel and perpendicular 
to Q, i.e. 

jP Q)Q 

P\\Q = Q2 ' (5) 

P±Q = P P\\Q. (6) 

We note that the QLBE described by Eqs. ([1]) and ([2]) has the structure of a translation-covariant master equation 
in Lindblad form, according to the general characterization given by Holevo |13l4l7| . This feature will be important 
below when applying the stochastic unraveling of the QLBE. 

B. Limiting forms 

Suitable limiting procedures reduce the QLBE to other well-known evolution equations, whose solutions are (at 
least partly) understood. These relations allow us to interpret the numerical solutions of the QLBE later on. At the 
same time, the stochastic simulation technique of the full QLBE permits us to study the range of validity of these 
approximate evolution equations. 



Classical linear Boltzmann equation 



To establish the connection to the classical linear Boltzmann equation one may consider the evolution of the diagonal 
elements w (P) = {P\p\P) in the momentum basis. As is shown in 0, [a, Q the incoherent part of the QLBE implies 
that 



dtw (P) = fdQ [Pr^ [P-Q^P)w{P-Q)- M'^ {P^ P + Q)w{P)] 



where the transition rates Af^' are given by 



M'^^P^P + Q) = / dk^\L{k^,P,Q)\ 
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Here a {Pf,Pi) = / {PfJPj) denotes the quantum mechanical scattering cross section. 

According to Refs. [3,[a[al, Eos. ([7]) and ^ agree with the coUisional part of the classical linear Boltzmann equation 
jl8| . In addition, it is argued in Q that the solution of the QLBE becomes asymptotically diagonal in the momentum 
basis for any initial state poj that is (Pje'^'Vol-F'' 7^ P) — >■ as i — >■ 00. It follows that the QLBE asymptotically 
approaches the classical linear Boltzmann equation for the population dynamics in momentum space. This fact will 
be important below when analyzing the diffusive behavior exhibited by the numerical solution of the QLBE. 

2. Pure colUsional decoherence 

The complexity of the QLBE reduces considerably if one assumes the test particle to be much heavier than the 
gas particles. By setting the mass ratio m/M equal to zero the Lindblad operators in ^ no longer depend on the 
momentum operator P of the tracer particle, so that the fc^-integration in ^ can be carried out |a,@]- The QLBE 
then turns into the master equation of pure coUisional decoherence 

MM, 
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where G (Q) denotes the normalized momentum transfer distribution and Feff is the collision rate of the gas environ- 
ment, defined by the thermal average 
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It leads to a localization in position space as can be seen by neglecting the Hamiltonian part in (jlOp for large M. The 
solution then takes the form 



{X\p{t)\X') = e-^(^-^'>(X|p(0)|X' 



(12) 



The decay rate of spatial coherences is given by the localization rate F (AX) > which is related to the momentum 
transfer distribution G (Q) by 



F{X- X') = Feff 



1- f dQGiQ)exp(^Q- {X - X' 
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The localization rate can be determined from the microscopic quantities as 



F (X - X') = Feff - 27rngas / dvfi{v)v dcosO \f {cos9; E^in) 
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where 9 denotes the scattering angle. Here we have assumed isotropic scattering, so that f {pf,Pi) = 
/(cos(pj,Pj) ;-Ekin =pf/2?Ti). Equation ([Tl)) will allow us below to predict the decoherence dynamics exhibited 
by the numerical solution of the QLBE in the limit M ^ m. 

3. Born approximation 

Another simplification results when the interaction potential V {x) is much weaker than the kinetic energy E = 
p'^ /2m. One may then replace the exact scattering amplitude / by its Born approximation fg, which is determined 
by the Fourier transform of the interaction potential, 

fB {Pf -P^) = -£kJ d-^ (-) exp (-^^^^f^) ■ (15) 

The approximated scattering amplitude therefore depends on the momentum transfer Pf^Pi only, so that the function 
/ in Eq. ([3]) is not operator-valued anymore. Taking ^ to be the Maxwell-Boltzmann distribution 

one may then perform the fe^-integration in ©, such that the dissipator V defined by Eq. ^ becomes [ll, y^lMSl 

VbP ^ J^Q (e'^'-^/'^is (P, Q) pL^B (P' Q) e"''^-^/'' 

-i{p,4(P,Q)LB(P,Q)}) . (17) 

Here the Lindblad operators contain the functions Lb (P, Q), given by the expression [l|, 0, Q 
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where as (Q) = |/b (Q)| denotes the differential cross section in Born approximation and /3 = 1/kT is the inverse 
temperature. The QLBE in Born approximation defined by Eqs. (|17p and (fT8)) was first proposed by Vacchini in 
Refs. [l|, |2[. As already mentioned, its solution may be obtained numerically by the stochastic simulation algorithm 
constructed in Ref. [12| . 

^. The limit of quantum Brownian motion 

The quantum Brownian motion or diffusion limit applies when the state of the testparticle is close to a thermal 
equilibrium state and when its mass is much greater than the mass of the gas particles [5|, |6| . The momentum transfer 
Q is then small compared to the momentum of the tracer particle. As discussed in [2l|, this permits the expansion of 
the Lindblad operators in (1^ up to second order in the position and momentum operators. This expansion yields the 
Caldeira-Leggett equation [ll|, [22] in the minimally extended form as required to ensure a Lindblad structure [y, [2l[ , 



^^ = ^[^-^] + l|[^'^^ + ^^]-T^ 



th 
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Here, A^y^ = 27r/i^/3/M gives the thermal de Broglie wave length, and 7 is the relaxation rate. It is remarkable that 
the derivation leads to a microscopic expression for the latter [6|, 
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The velocities of the gas particles are here assumed to be Maxwell-Boltzmann distributed, and the scattering to 
be isotropic so that the amplitude / depends only on the scattering angle 9 and the modulus of the momentum 
p= |pj I = |p f I • The integration variable u = p/pp denotes the momentum in dimensionless form, where pp ~ ^J2m/ (3 
is the most probable momentum at temperature T = 1/ {ksP). 

III. MONTE CARLO UNRAVELING 

To solve the QLBE we now employ the Monte Carlo wave function method |7Hll|. The underlying idea of this 
approach is to regard the wave function as a stochastic process in the Hilbert space of pure system states, with 
the property that the expectation value p{t) = ¥,[\'4> {t)){ijj {t) |] satisfies a given Lindblad master equation, dtp = 

(ifi.)~^[H, p\ + 'Yliii^iP^i ~ M'-I '-'' '^J')' ^'^y process with this property is called an unraveling of the master equation. 
An appropriate stochastic differential equation defining such a process is given by |ll| 



|dV*) - -|Heff|^t)dt + i^7.||U|V^t)f|^t)di 



u 

where Hoft represents the non-Hermitian operator 



Sl^a-i*'»^».('). (^') 



Hc« = H-'^Y.-'rllU^ (22) 

i 

The random Poisson increments ANi (t) in Eq. (|2ip satisfy the relations 

dN, (t) diVj (t) ^ S,,dN, it) , (23) 

and their expectation values are given by 

E[dN,it)] = 7jU|V't)|pdt. (24) 

The Monte Carlo method consists of generating an ensemble of realizations {IV'a (*))} of the process defined by the 
stochastic differential equation (PT|) , and of estimating the density matrix p (i) through an ensemble average |7Hll| . 

In the following, we briefly summarize a general algorithm which is often used for the numerical implementation 
of the stochastic differential equation (|2T|) [3]. This method forms the basis for the stochastic algorithm presented 
below, which extends the procedures presented in Ref. |12l |. 

A. The general algorithm 

We start with the normalized state \ipt) which has been reached through a quantum jump at time t (or is the initial 
state). Subsequently, the state follows a deterministic time evolution which is given by the nonlinear equation of 
motion 

dtii^t) = -'-H,s\iJt) + lY.^^\\^A)rH't) (25) 

i 

with the formal solution 

I , \ exp(-zHcffr/;t)|V;t) ,„„^ 

||exp(-zHoffT/ft)|V't)|| 

The probability for a jump to occur out of this state is characterized by the total jump rate 

r(^*) = ^E^[d^^W]=ETdlU|V*>lP. (27) 



It allows one to evaluate the corresponding waiting time distribution W {rlipt), the cumulative distribution function 
representing the probability that a jump occurs in the time interval [t, t + t], 

W{T\i^t) = l-!|exp(-zHeffT/fi)|7At)||^ (28) 

In practice, a realization r of the random waiting time can be obtained by the inversion method, i.e. by numerically 
solving the equation 

r, = |lexpHHeffT/?i)|V't)|P (29) 

for T, with ry a random number drawn uniformly from the interval [0, 1]. At time t + r a discontinuous quantum jump 
occurs, i.e. the wave function |V't+r) is replaced according to 

m+r) -^ TTT-r-, yT- (30) 

The corresponding jump operator, labeled by the index i, is drawn from the probability distribution given by the 
ratio of the jump rate Ti {ipt+r) = E [dNi {t + r)] /dt of the Poisson process Ni (i) and the total jump rate T {tpt+r), 

Prob(.|^.,.) = ^^.^4yf^. (31) 

B. Unraveling the QLBE 

We now adapt the Monte Carlo method to solve the QLBE, which is characterized by the family of Lindblad 
operators e^^'^'^'^L {kj_, P, Q). For this purpose, the index i is replaced by the continuous variables Q and k±, and 
the sums over i are substituted by the integrals 

V ^ [ dQ [ dfcx . (32) 

Although the procedure is straightforward, we repeat here the main steps since the obtained formulas are required 
for reference later on. 

The Monte Carlo unraveling of the QLBE is described by the stochastic Schrodinger equation 

\d^t) = ~l:^cffHt)dt+l f dQ f dk^\\L{k^,P,Q)\i^t)f\ijt)dt (33) 

^^i^. "'- [ ||L(fc„P,g)|^.)|i - 1^'*> j ^^^^'^^ ^*) ' 

where the effective Hamiltonian has the form 

Heff=H-^/dQ/ dfcxLt(fcx,P,Q)L(fcx,P,Q). (34) 

The Poisson increments in (|33|) have the expectation values 

E[dNQ,kAt)] = ||L(fcx,P,Q)K't)fdi, (35) 

and satisfy 

dNQ,k^ it) dNQ,^k'^ (t) = <5(3) ^g „ g/) ^(2) (^^ _ ^,^^ j^^_^^ (^) _ (gg) 

These relations represent the continuous counterpart of the discrete set of equations (j23p . The deterministic part of 
the Monte Carlo unraveling is generated by the nonlinear equation 

dtli^t) = -j^cffH't) + \ f dQ f dfcx||i(fcx,P,Q)|V^.)|PlV't), (37) 

" ^ Jr3 Jq± 



whose formal solution is given by Eq. ([26]). The jump probability is determined by the rate 

rCV'O = [ dQ [ dk^\\L{kj_,P,Q)\^t)\\\ (38) 

and a realization of the random waiting time r is obtained by solving Eq. (j29p for r with the effective Hamiltonian 
([M]) . The jump at time i + r is effected by 

\\L(kj_,P,Q)\ipt+T)\\ 
where the continuous parameters fcx and Q characterizing the jump operator arc drawn from the probability density 

Pvoh{k^,Q\^JJt+r) = —y^-—A\Lik^,P,Q)\^t+r)f. (40) 

C. Unraveling the QLBE in the momentum basis 

The implementation of the above algorithm is particularly simple when the initial state is a discrete superposition 
of a finite number of momentum eigenstates [12| , 

JV JV 

1^ (0)) = ^ a, (0) \P, (0)) , with J2 l«^ (O)l' = 1 • (41) 

4=1 1=1 

Due to the translation-covariance of the QLBE the Lindblad operators have the structure ef'^'^I^L {k±, P, Q). This 
implies that the effective Hamiltonian is a function of the momentum operator only, so that the deterministic evolution 
of (|4T|) affects solely the amplitudes of the superposition, that is 

N 

|V'(i))=^a,(t)|P,(0)). (42) 

4=1 

The jumps, on the other hand, cause a translation of the momentum eigenstates and a redistribution of the amplitudes, 

N 

,^Q->^in^ (fcx, P, Q) \i, (t)) =Y.< (^) \P- + Q) • (43) 

This shows that the quantum trajectory |'0(i)) remains a superposition of -/V momentum eigenstates at all times. The 
stochastic process therefore reduces to a process in the finite-dimensional space of the amplitudes on and momenta 
Pi. Here the momentum eigenstates are taken to be normalized with respect to a large volume $7, such that they 
form a discrete basis, {Pi\Pj) = Sij. 

In the following it is convenient to work with dimcnsionless variables 

C/^^, K^-^, W^^^, (44) 

Mvp in^vp mvp 



where the scale is given by the most probable velocity of the gas particles vp = ^J^kj^Tjrn. Note that W^j_, being 
proportional to fcj^, lies in the plane perpendicular to K. The quantum trajectories are then represented as 



JV N 



\^ (i)) = Y. "■' (^) 1^* (*)) ' ^ith E I"' (*)!' = 1 ■ (45) 
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Before discussing the unraveling of the QLBE in more detail, let us evaluate the jump rate psp for momentum 
eigenstates, j-f/'t) = |P). This quantity appears frequently in the algorithm described below. By inserting jf/^t) = \P) 
into Eq. psp . one obtains 

r(P) = / dQ / dfci|L(fci,P,Q)|2. (46) 

JR3 jQi 



Noting Eq. ([8|), one finds that the jump rate agrees with the total coUision rate for a particle with momentum P, 

r (P) = f dQAr' {P^P + Q) . (47) 

It follows that r (P) = r (P) is a function of the modulus of P only, since the collision rate must be independent 
of the orientation of P for a homogeneous background gas. Upon using the dimensionless quantities (j44p . and after 
inserting (j3]) for L, as well as the Maxwell-Boltzmann distribution (fT6|) . one finds 



r(;7) == f dK f dW^giW^,U,K)p,^{K)p„,,{WA 



(48) 



with 



g{W^,U,K) = 
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The densities pa-^^ (K) and Pay^- (W ±) denote three- and two-dimensional normal distributions, respectively. 






1 



1 

2'K(7\Y 



3/2 

exp 



exp 



2^ 



2^ 



(50) 



with variances ax = v2 and aw = l/v2- 

The integral (|48|) can be evaluated numerically using a Monte Carlo method with importance sampling [23] ■ For 
this purpose, one draws n samples Ki from the normal distribution p„j^ {K) and computes orthonormal vectors en 
and e2i which arc orthogonal to Ki, i.e. en ■ Ki = 0, 62; • Ki = 0, en ■ e^i = using the Gram-Schmidt method. As 
a next step, n further samples (ui, v^i) are drawn from the two-dimensional Gaussian distribution PaT,^-, which yields a 
sample of scaled momentum vectors W i^ — UiCn + Vie2i- The jump rate (|48p is then approximated by the average 
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V{U) ~ -yg{W^^,U,K,) . 
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Let us now discuss in more detail the unraveling of the QLBE in the momentum basis. To this end, suppose the 
state 
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\^{t)) = Y.^dt)\udt)) 



(52) 
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was obtained through a quantum jump at time t. As mentioned above, the effective Hamiltonian p4p depends on the 
momentum operator only, so that the momenta Ui stay constant during the deterministic evolution. The propagation 
of the state ([5^ with the non-Hermitian operator ([M)) thus yields [l2| 
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Here the weights have the form 



a, [t + r) = -1 exp (-l^MvpUjA exp (-^F (Ui) ] a, (t) , 



with the normalization 
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As a next step, one must evaluate the waiting times r. For this purpose, consider the expression 
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By using the definition of Hcg p4p. the fact that the two summands in Hofr commute, and the jump rate (j48p . this 
yields 
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exp(-zHeffr/^)lV't)f = ^ |a, (i)|' exp (-rF (C/,)) 
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It follows from (|29p that samples of the waiting times t are obtained by numerically inverting the non-algebraic 
equation 



N 



= ^|a,(i)|'exp(-Tr((7,)), 



(58) 



with rj drawn from the uniform distribution on [0, 1]. 

To be able to carry out the quantum jumps, we have to determine the momentum parameters K and VFx, which 
characterize the jump operator. These vectors are obtained by sampling from the probability distribution (j40p . Upon 
inserting states of the form ([53| , Eq. (|40|) becomes 
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This distribution is a mixture of the probabilities 
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and the probability densities 
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In order to draw a sample from the mixture (|60p . one may proceed as follows [12[. First, an index i is drawn from the 
probabilities (|B(I|. Then, the momenta if and VFj^ are drawn from the probability distribution Proh {W±,K\Ui) 
using a stochastic sampling method, such as the Metropolis- Hastings algorithm [23[. 

Having the momenta W j_ and K at hand, one can now perform the quantum jump. According to Eq. ()39p . the 
state ([53]) is transformed as 

\i>{t + T)) ^ AA-^expf ^mi;^l<:-X]L(VFi,U,K)^a,(t + r)|[7,(t)) 
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where the normalization A/" is determined by 
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FIG. 1; CPU time of the Monte Carlo unraveling as a function of the number A'' of basis states involved in the initial 
superposition. The curve is almost a straight line with slope a ~ 1.1 in the logarithmic plot. It follows that the CPU time 
grows almost linearly with A*', that is T oc A^^'^. 



This shows that the momentum eigenstates are shifted 

while the weights are redistributed as 

ttj {t + r) -> a[ it + r) = Xitti (i + r) , 
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where the factors Xi arc given by x^ = A/" ^L {W±, Ui, K). Upon using the explicit form ([3]) of L, and by inserting 
the Maxwcll-Boltzmann distribution (I16p . we find 
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According to Eq. (|64p . the momentum eigenstates are all shifted with the same momentum K in a. quantum 
jump. This fact is decisive for the numerical performance of the algorithm, since it implies that the time consuming 
Metropolis-Hastings algorithm must be applied only once for alH G {1, . . . , N}. This suggests that the algorithm can 
be applied also to initial states which are superpositions of many momentum eigenstates. 

This fact is substantiated by the numerical analysis depicted in the logarithmic plot of Fig. [T] Here, the CPU time 
of the above algorithm is shown as a function of the number N of basis states involved in the initial superposition 
state 1-0 (0)) = X]i=i \Ui) /\fN , with Ui — (0, 0, i). The simulation is based on 10^ quantum trajectories in each run. 
The curve shown in Fig. [1] is almost a straight line with a slope a ~ 1.1, implying that the CPU time T grows almost 
linearly with TV, T oc N^-^. 

We conclude that the Monte Carlo unraveling can be implemented for initial superposition states that are composed 
of a large number of momentum eigenstates (say, on the order of 10^ to 10'^). This implies that one may choose even 
well localized initial states and consider scenarios where a particle crosses a slit or a grid. The following sections 
present numerical results obtained with such kinds of states. 



IV. SIMULATION RESULTS 



We proceed to apply the stochastic algorithm to two different types of scattering interactions with the surrounding 
jas particles. Specifically, we consider the simplest possible scattering process (s-wavc hard-sphere scattering) as well 
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as the case of a general potential, which is treated exactly through partial wave decomposition. Having discussed 
the determination of the scattering amplitudes, we start out with the simulation of short-time effects. At first, the 
loss of coherence of an initial superposition of two momentum eigenstates is measured, followed by the treatment of 
superpositions of spatially localized wave packets. The latter permits in particular to extract the localization rate 
discussed in Sect. IIIB 21 As a further example of a decoherence process, counter-propagating localized initial states 
are considered which lead to the formation of interference patterns. In the course of the evolution, fringe visibility is 
lost, so that the interplay between coherence and decoherence can be demonstrated. 

We then discuss long-time effects which exhibit a classical counterpart, starting with energy relaxation and the 
approach to thermal equilibrium. Then, the spread in position of initially spatially localized states is measured, 
allowing us to observe a transition from quantum dispersion to classical diffusion. 

As discussed in Sect. IIIBl the QLBE has several limiting forms for some of which analytical solutions are known. 
This permits to demonstrate the validity of the numerical results and to verify the limiting procedures discussed in 
[6|. Further simulations correspond to situations where the full QLBE is required. This way physical regimes are 
entered which have not been accessible so far, such as decoherence phenomena where the mass of the test particle is 
comparable to the mass of the gas particles. 

A. Scattering amplitudes 

1. S-wave hard-sphere scattering 

In s-wave hard-sphere scattering the particles arc assumed to be hard spheres with radius R, and the kinetic energy 
to be sufficiently small, pR <C h, such that only the lowest partial wave contributes. In this case the scattering 
amplitude is independent of the scattering angle and the kinetic energy, |/ (cos0; _Ekin)| = R^- For a constant 
cross section one can do the fej^-integration in the QLBE ©. The equation then coincides with the QLBE in Born 
approximation ((T7)) . such that the numerical results with this interaction should agree with the stochastic algorithm 
of Breuer and Vacchini [i2[. In the s-wave examples presented below the system of units is defined by setting h = 1, 
M = 1 and i? = 1; the temperature is chosen to be fc^T = 1 and the gas density is set to one, rigas = 1- 

An important ingredient for implementing the Monte Carlo unraveling is the jump rate F iU) presented in Eq. (j48p . 
It is obtained numerically by Monte Carlo integration with importance sampling (|5ip . based on n = 10^ steps. We 
find that the collision rate grows linearly for large momenta, while it saturates for vanishing [/ at a value close to 

To = rigas w/347riT!^ (67) 

in agreement with the analytical prediction in |12| . 

2. Gaussian interaction potential 

Generic scattering processes are characterized by many partial waves with energy dependent scattering phases. 
To illustrate the treatment of this general scattering situation, we choose the scattering amplitudes defined by an 
attractive Gaussian interaction potential 



V{r) = -FoCxp(^-^j . (68) 

The corresponding scattering amplitude in Born approximation is obtained from (jlSp . which yields 

/b(cos6I,p) = J- exp — -^[l-cos6i] . (69) 



ft2 -^y ^2 

While this approximation is reliable only for weak interaction potentials, Vq ^ £'kin, the exact scattering amplitudes 
require the energy dependent partial scattering amplitudes /; [2J], 

oo 

/ (p, cos e) - ^ (2; + 1) /, (p) Pi (cos e) , (70) 

with Pi the Legendre-polynomials. 
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FIG. 2: The first four pliase sliifts 5i as a function of tlie relative momentuni p for the interaction potential (|68|l and mass ratio 
m/M = 1. Their asymptotic behavior is in agreement with (|71[) . (There are two bound states for I = 0,1, one for I = 2, and 
no bound state for / = 3.) 

The /; (p) = (fi/p) exp{i5i) siii(Si) are related to the partial wave phase shifts Si, which can be computed numerically 
by means of the Johnson algorithm [25| for a given interaction potential. If the kinetic energy is large compared to 
the potential, V <C p^ /2M, the partial waves are hardly affected by the collision, so that the scattering amplitudes 
and phases vanish, Si {p -^ oo) = 0. For small energies, on the other hand, they behave as [2J] 



Slip) 



71; TT - aip 



„2i+l 



for p^ , 



(71) 



with ai the scattering lengths and ni G No. According to the Levinson theorem [24| . the integer ni equals the number 
of bound states with angular momentum I. Figure [2] shows the first four phase shifts for Vq = 20, d = 1, and fi, = 1, 
in agreement with the Levinson theorem (|71|) and with the expected high energy limit. 

For the simulations presented below it is sufficient to include the first 30 partial waves when evaluating the scattering 
amplitudes ((70|) . In particular, this ensures that the optical theorem is satisfied |24| . 

Figure [3] shows the numerically evaluated jump rate F (U) for the Gaussian interaction potential. It is obtained by 
a Monte Carlo integration of Eq. (|48| with importance sampling with n = 10^ steps. The jump rate is given in units 
of the collision rate as defined by the thermal average PT|) . The simulation shown by the solid line in Fig. [3] is based 
on the exact scattering amplitude, while the dashed line corresponds to its Born approximation. One observes that 
the two results differ drastically, in particular for large interaction potentials Vq, while they tend to agree for large 
momenta p, where the Born approximation is more reasonable. 

The Gaussian interaction potential is applied in several examples below. In these cases, the system of units is 
defined by setting ?i = 1, m = 1 and rf = 1; moreover, we chose ksT = 1 for the temperature of the gas environment 
and the gas density is set to unity, ngas = 1. 



B. Decoherence in momentum space 

We now apply the Monte Carlo algorithm to the analysis of decoherence effects in momentum space. For this 
purpose the initial state is taken to be a superposition of two momentum eigenstates, 



IV-W) = «(0)|C/(0))+/3(0)|V(0)), 



(72) 
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FIG. 3: Jump rate F as a function of the momentum U assuming the Gaussian interaction potential with Vb = 1 (left) 
and Vo = 20 (right). The sohd hue corresponds to the exact scattering amplitude, while the dashed line gives the Born 
approximation. As one expects, the results deviate strongly for large interaction potentials. 
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FIG. 4: Semi-logarithmic plot of the "coherence" C (t) defined in (|75|l for the state ((72} with Uo = VG- The interaction is 
described by the Gaussian potential (|68|) with Vb = 1 (left) and Vb = 20 (right). The solid line is obtained using the exact 
scattering amplitude, while the dashed line corresponds to the Born approximation. The predictions of the decoherence rates 
differ substantially in case of the large interaction potential with Vb = 20. 



which are assumed to have the form U (0) ~ ~V (0) = {Uq, 0, 0). 

Since the states \U (0)) and | V (0)) are genuine momentum eigenstates, any collision necessarily leads to an orthog- 
onal state. It follows that the coherences are expected to decay exponentially 



\{UiO)\pit)\ViO))\ 
|(C/(0)|p(0)|V(0))| 



-r{Uo}t 



(73) 



with the decay rate given by the total collision rate F (Uo). 

Alternatively, one may view the states \U (0)), \V (0)) as representing states which are well localized in momentum 
space, but with a finite width greater than the typical momentum transfer. Here a suitable measure for the degree 
of coherence is the ensemble average of the coherences exhibited by the individual quantum trajectories \ip (t)) |12| . 
that is 



C{t) = E 



\{u{t)\^it)){^it)\vm 



|(c/(o)|p(o)|y(o))| 



(74) 



To evaluate this term, recall that the quantum trajectories remain in a superposition of two momentum eigenstates 
all the time, so that \ip (t)) has the form \tjj (t)) = a (t) \U (t)) + /3 (t) \V (t)). By inserting this expression into equation 
dZl), one finds ^ 



C{t) = 2E[|a(t)/3*(t)|] 



(75) 
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Figure |4] shows a semi- logarithmic plot of the "coherence" C (t) for the Gaussian interaction potential, choosing an 
initial momentum Uq = y/6, equal amplitudes a (0) = /3 (0) = l/"\/2, and the mass ratio M/m = 1. The left-hand side 
represents a weak interaction potential, and the right-hand side a strong one. In the latter case, the result obtained 
with the exact scattering amplitude (solid line) differs markedly from the corresponding Born approximation (dashed 
line). The simulation is based on 5 x 10'^ trajectories. 

This result shows that the full QLBE ([2]) may lead to physical predictions which deviate significantly from the 
ones obtained with the QLBE in Born approximation (J17p if the interaction potential is sufficiently strong. A similar 
conclusion is drawn below, when studying relaxation rates. 

The design of experimental tests for decoherence effects in momentum space is a challenging task 0, [2^, [23] ■ Such 
a setup would have to provide a source of states with momentum coherences (as in non-stationary beams), and it 
would require an interferometric measurement apparatus able to detect these coherences. A further difficulty lies 
in the inevitable presence and dominance of position decoherence. During the free evolution a superposition state 
characterized by two different momentum values will evolve into a superposition of spatially separated wave packets, 
which is affected by decoherence mechanisms in position space |6| . 

Position decoherence, in contrast, has already been observed experimentally in fuUerene interference experiments 
[28| . The following section therefore focuses on the prediction of spatial decoherence effects based on the Monte Carlo 
unraveling of the QLBE. 

C. Decoherence in position space 

1. Measuring spatial coherences 

In order to quantify the loss of spatial coherences, i.e. the off-diagonal elements in position representation, 
p(X,X'] = {X\p\X') one must assess p(X,X'] given the quantum trajectories in the momentum representa- 
tion, \ip (t)) — X^i^i '^j W \^j W)- For this purpose, it is convenient to express the position variable X in units of 
the thermal wavelength Ath = \/2'nh? /mkBT , 

S ^ ^. (76) 

The spatial coherences are then obtained by taking the ensemble average of the coherences of the individual quantum 
trajectories, that is 

p{S,S',t) = E[{S\i:{t)){^{t)\S')] . (77) 

By inserting the momentum representation of |i/' (i)) into this expression, we find 

p{S,S',t) (78) 

— 3 ^ E ^a, (t) al (i) exp f ^Af ^^Ath [S ■ U, (t) -S'-Uk (i)] j 



(2^)^ 



j.k 



which allows us to compute the time evolution of the coherences (|77p by means of the amplitudes aj (t) and the scaled 
momenta Uj (t). 

A typical application might describe a particle passing through an interferometer, where it is spatially localized in 
one spatial direction, and is characterized by an incoherent distribution of momenta in the other two directions. From 
now on, we therefore restrict the discussion to initial states of the form 

N 

1^ (0)) = E "^- (0) \U^ (0) ' ^ (0) ' ^ (0)) ' (79) 

where \Uj (0) ,V (0) ,W (0)) denote scaled eigenstates of the momentum operator P = (P^,, P^, P^). By taking N to be 
sufficiently large, Eq. (j79p may represent states which are localized in one spatial direction. Due to the conservation 
of momentum superpositions, the ensuing quantum trajectories have the structure 

JV 

\i,{t)} = Y.a,{t)\U,{t),V{t),W{t)). (80) 
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FIG. 5: Evolution of the density matrix in the position representation for an initial superposition of two Gaussian wave packets, 
obtained by solving the three-dimensional QLBE for s-wave hard-sphere scattering. The spatial coherences p{x/A.th,y/f^th) 
are expressed in units of the thermal wavelength. 



The assessment of spatial coherences (|78p can be simplified in this case by focusing on the coherences in x-direction, 

p([5,0,0],[5',0,0],i) (81) 



N 






a, [t) al (i) exp ( -Mvp Ath [SUj (i) - S'Uk (t)] 



To visualize the evolution of the density matrix in position representation, we consider an initial superposition of 
two resting Gaussian wave packets, with scaled mean positions (S)i^2 = ±1-2 and width cri.2 = 0.2 (in units of Ath). 
This state may be written in the form (|80|) by using a finite-dimensional representation of the corresponding Fourier 
transform. Figure [S] depicts the ensuing evolution of the matrix elements (|8ip , obtained by solving the QLBE under 
the assumption of s-wave hard-sphere scattering and equal masses m = M . It shows four snapshots of the density 
matrix for the scaled times iFo G {0,1/3,2/3,4/3}. The simulation is based on 10'^ realizations of the stochastic 
process and the state is represented using N = 55 momentum eigenstates. 



2. Measuring the localization rate 



As discussed in Sect. IIIB2[ the QLBE simplifies to the master equation of pure collisional decoherence if one 
assumes the tracer particle to be much heavier than the gas particles. In this model the decay rate F of spatial 
coherences is a function of the distance x = \X — X \ only; it does not depend on the particular matrix elements of 
the state, see Eq. (fT4|) . Hence, the decoherence process is completely characterized by the localization rate F (x). 

By evaluating the decoherence rates for various mass ratios and initial states, we found that this behavior is observed 
even in regimes where the QLBE does not reduce to the master equation of collisional decoherence. This suggests 
that the decoherence dynamics of the QLBE is generally characterized by a one-dimensional function F {x). 

Figure [S] shows the localization rate for the Gaussian interaction potential with Vq = 1 and the mass ratios 
Mjm — 100 (left) and M/m = 1 (right). The dots give the decay rate as evaluated from ([ST|) . obtained by 5 x 10"' 
realizations of the Monte Carlo unraveling of the QLBE. The solid line represents the localization rate of collisional 
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FIG. 6; Decay rate of the spatial coherences as a function of the wave packet separation for a Gaussian interaction potential 
and the mass ratios M/m — 100 (left) and m — M (right). The solid line shows the prediction of pure coUisional decoherence, 
Eq. (|14p . and the dots give the result of the stochastic simulation of the QLBE. One observes that the predictions of the two 
models agree for M ^ m, while they deviate for m — M. The localization rate saturates in all cases at the average collision 
rate Fefi. We note that the decay rate of the QLBE does not vanish for m = M as x -^ Q, since there is a loss of the populations 
due to diffusion. 

decoherence (J14p . calculated by numerical integration. As expected, one finds an excellent agreement between the 
predictions of coUisional decoherence and the solution of the QLBE if the test particle mass is much larger than the 
gas mass, M/m = 100. 

Moreover, it turns out that the results of the two models do not differ substantially even for equal masses m = M . 
This holds in particular for large distances, where the decay rates converge to the average collision rate Tcs (in all 
cases). Indeed, in this limit one collision should be sufficient to reveal the full 'which path' information, so that a 
saturation at T^s is expected. For equal masses the prediction of the QLBE does not tend to zero in the limit of small 
distances, F (0) > 0. This is due to the contribution of quantum diffusion, which is more pronounced when the test 
particle is lighter. 

D. Interference and decoherence 

To illustrate the interplay between coherent and incoherent dynamics, let us study how the formation of interference 
patterns is affected by the interaction with the background gas. To this end, consider the scenario depicted in 
Fig. [T] Here the x-component of the three-dimensional initial state is prepared in a superposition of two counter- 
propagating minimum-uncertainty wave packets ^1,2, while the other two components have a definite momentum. 
The wave packets start overlapping in the course of the evolution, and their interference leads to oscillations of the 
spatial probability density p{x,x,t), with a period given by the de Broglie wavelength AdB associated to the relative 
momentum between the minimum-uncertainty wave packets. Besides this coherent effect, one observes an increasing 
signature of decoherence, the gradual loss of fringe visibility; this becomes evident in particular in the bottom panel 
of Fig. El 

Figure [7| is obtained by the Monte Carlo unraveling of the QLBE, assuming s-wave hard-sphere scattering and a 
mass ratio M/m = 100. The parameters of the simulation are conveniently expressed in terms of the de Broglie 
wavelength AdB and the scattering rate Fq (|67p . which serve to define the dimensionless variables 



5dB 



X 

AdB 



f/dB = 



P 



M AdBFo 



(82) 



In this system of units the position and momentum expectation values of the coherent states ^12 read as (SdB)i,2 = 
=F15 and (UdB)i.2 = ±0.9; their width is characterized by the standard deviation cri,2/AdB ~ 4, and the de Broglie 
wavelength is fixed by setting AdB/Ath — 2.5 x 10^^. The figure shows three snapshots of the populations of the 
density matrix for the scaled times Fg (ioi ^1,^2) ~ (0, 9, 18). The simulation is based on 2.5 x 10^ realizations of the 
stochastic process. 

As mentioned above, an important quantity to characterize the loss of quantum coherence is the fringe visibility, 
which we define here as the difference between the central maximum and the neighboring minimum divided by their 
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FIG. 7: Evolution of the position diagonal elements of the density matrix p {x/Xds, s/Ads) for an initial superposition of two 
counter-propagating minimum-uncertainty wave packets. The figure, obtained by solving the QLBE for s-wave hard-sphere 
scattering, shows three snapshots of the dynamics at times To (^0,^1,^2) = (0, 9, 18). The increasing influence of decoherence 
manifests itself as a reduced visibility of the interference fringes as time progresses. 

sum. In the last snapshot, shown at the bottom of Fig. [71 one extracts a visibility of V = 55%. To understand this 
result quantitatively, let us estimate the decay rate of the visibility by means of the integrated localization rate. 



V(t) = cxp(- /" dr'F[S'(r')]T'j V(0) , 



(83) 



where S (r) ~ \ (X (t)) 1 — (X (t))2 | /Ath denotes the distance of the minimum-uncertainty wave packets in units of the 
thermal wavelength at time r. By noting that the wave packets move in absence of an external potential, one finds 



Sir) = ^(l(SdB)i 
Ath 



3dB;2 



r|(U 



dB/l 



(U 



dB;2 



I) 



(84) 



Since the tracer particle is much heavier than the gas molecules the dynamics described by the QLBE should be well 
approximated by the master equation (fTO|) of pure collisional decoherence. In this case F is described by the formula 
([M]) . which can be evaluated analytically in the case of s-wave hard-sphere scattering. 



F{S) = 2V^ngasi?\;3[4-S'^^exp(-47r52)erfi(2V^S')] , 



(85) 



where erfi(a;) = — ierf (ix) denotes the imaginary error function. 

A prediction for the visibility (|83p may then be obtained by a simple numerical integration. This yields V (^2) — 56%, 
in good agreement with the value of V = 55% obtained by the stochastic solution of the full QLBE. 



E. Relaxation and thermalization 



We now study the long-time behavior of the energy expectation value. As discussed in [6| any solution of the 
QLBE will approach the canonical thermal state asymptotically. The kinetic energy in the simulation must therefore 
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FIG. 8: Energy relaxation as obtained by using the exact scattering amplitude of the Gaussian interaction potential, compared 
to the corresponding Born approximation and the solution of the CL equation (|20|l . We choose the potential strengths Vo = 1 
(left) and Vo = 20 (right) and the mass ratios M/m = 1 (top) and Mjra — 10 (bottom). The correct equilibrium values (|86p 
of 3/2 (top) and 3/20 (bottom) are obtained in all cases. One observes that the exact result agrees with the solution of the 
CL equation for a heavy tracer particle (bottom), while it deviates strongly for M = m (top). The Born approximation gives 
reliable results only when the kinetic energy is much greater than the potential one (top left), which indicates that it may 
underestimate the energy relaxation even for weak interactions if the test particle has a large mass (bottom left). 



converge to the thermal energy 3/ (2/3). Expressed in dimensionless units, this means that 



(86) 



with 7 the relaxation rate. 

If the state is close to thermal, and if the tracer particle is much heavier than the gas particles the QLBE reduces to 
the Caldeira-Leggett (CL) equation in Lindblad form ()19|) . The corresponding energy behavior is then well understood 

@,[ni. 



(U^), = (U2)eq+((U2),„-(U2)cq) 



,-47t 



Moreover, the relaxation rate 7 can be expressed in terms of the microscopic quantities, see Eq. ((201 
can be evaluated analytically for the case of a constant cross section [0] , 



(87) 
The integral 
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Figure [8] shows the energy relaxation exhibited by the stochastic solution of the QLBE for a weak and a strong 
Gaussian interaction potential, with Vb = 1 (left) and Vq = 20 (right) respectively. The solid line depicts the solution 
of the QLBE based on the exact scattering amplitude, while the corresponding Born approximation is represented by 
the dashed line; both simulations are based on 5 x 10^ trajectories. The initial state is here a momentum eigenstate 
with dimensionless eigenvalue Uq = VG (top) and Uq = VO.G (bottom), corresponding to mass ratios of M/m = 1 
and M/m =10, respectively. In case of a relatively large tracer mass (bottom) one obtains a good agreement with 
the prediction of the CL equation (|87p (open dots) . Here the relaxation rate was obtained by numerical integration 
of the right-hand side of Eq. ([^0]) . For equal masses, on the other hand, the results deviate noticeably (top). As 
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expected; all of the solutions converge to the correct equilibrium values, given by the scaled energies 3/2 (top) and 
3/20 (bottom). The Born approximation yields reliable results only in the situation depicted by the top left panel, 
where the kinetic energy is much larger than the potential. 

Again, we are led to conclude that the full QLBE ^ may give rise to predictions which deviate significantly from the 
ones obtained with the QLBE in Born approximation (|17p. This holds in particular for strong interaction potentials, 
where the corresponding scattering amplitudes are different. Furthermore, this result verifies that the expression ((20|) 
obtained in [5|, l6| yields the correct relaxation rate in the quantum Brownian limit. 

F. Diffusion 

As a final aspect we study the quantum diffusion process described by the QLBE. To this end, a localized initial 
state is prepared and the growth of the position variance is measured. Before discussing the numerical result, we 
summarize analytical predictions based on [6| . 

On short time scales, where the number of collisions is small, one expects the variance growth to be dominated by 
quantum dispersion. This implies that the variance growth is parabolic; for an initial state of minimum uncertainty 
one expects 

At time scales after which many collisions have occurred the variance growth is expected to be dominated by 
classical diffusion. The corresponding diffusion constant can be estimated by considering that the QLBE approaches 
asymptotically the classical linear Boltzmann equation asymptotically. The latter can be simplified, by taking the 
Brownian limit of heavy tracer particles with a momentum P close to the typical thermal value P/j ~ ■\/2Af//3 
[a, [23,1301 ■ Under these conditions, the classical linear Boltzmann equation ([7]) reduces to the Kramers equation [yj 

d,w{P) = .,^(^_[P,^(P)] + __^.(P)j (90) 

for the momentum distribution wJP), with rj the friction coefficient. The latter can be expressed in terms of the 
microscopic details of the gas [g, |31|, yielding ?; = 2j, with 7 the relaxation rate appearing in the Caldeira-Leggett 
equation, see Sect. Ill B4l Eq. ((201) . 

The Kramers equation predicts normal diffusion, i.e. a linear growth of the variance (t|- (t) ~ (j\ (0) + 2_Dt, with 
diffusion constant D = rjM/ j5 [32[- This leads to the prediction that 

^l{t) = ^x{^) + ^t (91) 

whenever the Brownian limit of the QLBE is applicable. For the case of a constant cross section, where 7 can 
be evaluated analytically (see Eq. (|55)) ). Eq. (|?T|) provides an analytical prediction for the diffusion constant. It is 
expected to be valid when M ^ m. 

The solid line in Fig. [9] shows the variance growth of the spatial populations, obtained by solving the QLBE for 
s-wave hard-sphere scattering and mass ratios M/m = 100 (left) and M/m = 1 (right). This stochastic simulation is 
based on 4 x 10"^ trajectories. The initial state is chosen to be a Gaussian with width a^ (0) /A^j^ — 1.6 x 10'"^ (left) 
and cr^ (0) /A2j^ = 1.6 x 10"! (right). 

Let us first focus on the left-hand side panel which corresponds to a very massive particle. Here the solution of 
the QLBE starts with a quadratic dependence at small times. The curvature is unrelated to that of free quantum 
dispersion (dashed line), Eq. (|89|. which is clearly due to the large number of collisions occurring on the time scale 
of the wave packet broadening. (Time is given in terms of the average period between collisions in our dimensionless 
units.) 

After a time corresponding to about 200 collisions the curve displays the straight line behavior expected for classical 
diffusion. A linear fit to this straight part, indicated by the dotted line, has a slope of approximately 7.5 x 10^^. This 
differs by about 12% from the analytical considerations presented above, where one expects a straight line of the form 
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FIG. 9: The solid lines give the time dependence of the spatial variance, as obtained by solving the QLBE for s-wave hard- 
sphere scattering and mass ratios of M/m = 100 (left) and M/m — 1 (right). Left: after an initial quadratic increase the 
variance displays the straight line behavior expected of classical diffusion. The fit displayed by the dotted line has a slope close 
to the diffusion constant predicted by the Kramers equation (relative error: 12%). For comparison, the dashed line gives the 
dispersive broadening of the initial wave packet in absence of a gas. Right: for a test particle with a mass equal to that of the 
gas particles the variance growth is dominated by classical diffusion even on the time scale of a single collision. 

For equal masses (right panel) one obtains a straight line already starting from small times, which indicates that 
classical diffusion dominates over quantum dispersion. The slope of about 3.2 x 10~^ implies that the diffusion constant 
is much greater for light test particles. However, these results cannot be related to Kramers equation, since the latter 
is valid only in the Brownian limit of large tracer masses. 



V. CONCLUSIONS 



We presented a stochastic algorithm for solving the full quantum linear Boltzmann equation given an arbitrary 
interaction potential. By exploiting the translational invariance of the QLBE it allows one to efficiently propagate 
superpositions of momentum eigenstates without increasing the dimension of the state space. Since the computation 
time scales almost linearly with the number of basis states, arbitrary states can be represented in practice, in particular 
spatially localized ones. This enables us to simulate many important physical processes, ranging from short-time 
effects, such as the loss of fringe visibility in interference experiments, to long-time relaxation and thermalization 
phenomena. 

For the cases of s-wave hard-sphere scattering and a Gaussian interaction potential, we analyzed the range of 
validity of different limiting forms of the QLBE, including the collisional decoherence model, the quantum Brownian 
limit and the classical linear Boltzmann equation. Moreover, we compared the solutions of the full QLBE to those of 
the simplified equation in Born approximation. Here it is found, for the above interactions, that the full QLBE may 
lead to physical predictions which deviate significantly from the ones obtained with the QLBE in Born approximation 
if the interaction potential is sufficiently strong. 

This method will find applications, e.g. in describing interference experiments with species, whose mass is smaller 
than or comparable to the mass of the gas particles. The existing methods are not able to quantify the loss of 
coherence in such a situation. Moreover, future studies might consider extensions of the discussed method to the 
recently developed quantum master equation for the collisional dynamics of particles with internal degrees of freedom 
p33l - l35| . Even though this equation is more involved than the QLBE, it is also translational invariant. Since this 
property is the main prerequisite for the present algorithm, it should be extensible to the quantum master equation 

of [si. 
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